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-2009 Kevin Rosenberg
  11 **
  12 **  This program is free software; you can redistribute it and/or modify
  13 **  it under the terms of the GNU General Public License (version 2) as
  14 **  published by the Free Software Foundation.
  15 **
  16 **  This program is distributed in the hope that it will be useful,
  17 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 **  GNU General Public License for more details.
  20 **
  21 **  You should have received a copy of the GNU General Public License
  22 **  along with this program; if not, write to the Free Software
  23 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24 ******************************************************************************/
  25
  26 #include "ct.h"
  27
  28
  29
  30 void
  31 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)
  32 {
  33   ImageFileArray vReal = im.getArray();
  34   ImageFileArray vImag = im.getImaginaryArray();
  35   unsigned int ix, iy;
  36   unsigned int nx = im.nx();
  37   unsigned int ny = im.ny();
  38
  39   // shuffle each column
  40   for (ix = 0; ix < nx; ix++) {
  41     Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
  42     if (im.isComplex())
  43       Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);
  44   }
  45
  46   // shuffle each row
  47   float* pRow = new float [nx];
  48   for (iy = 0; iy < ny; iy++) {
  49     for (ix = 0; ix < nx; ix++)
  50       pRow[ix] = vReal[ix][iy];
  51     Fourier::shuffleFourierToNaturalOrder (pRow, nx);
  52     for (ix = 0; ix < nx; ix++)
  53       vReal[ix][iy] = pRow[ix];
  54     if (im.isComplex()) {
  55       for (ix = 0; ix < nx; ix++)
  56         pRow[ix] = vImag[ix][iy];
  57       Fourier::shuffleFourierToNaturalOrder (pRow, nx);
  58       for (ix = 0; ix < nx; ix++)
  59         vImag[ix][iy] = pRow[ix];
  60     }
  61   }
  62   delete pRow;
  63 }
  64
  65 void
  66 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
  67 {
  68   ImageFileArray vReal = im.getArray();
  69   ImageFileArray vImag = im.getImaginaryArray();
  70   unsigned int ix, iy;
  71   unsigned int nx = im.nx();
  72   unsigned int ny = im.ny();
  73
  74   // shuffle each x column
  75   for (ix = 0; ix < nx; ix++) {
  76     Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
  77     if (im.isComplex())
  78       Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);
  79   }
  80
  81   // shuffle each y row
  82   float* pRow = new float [nx];
  83   for (iy = 0; iy < ny; iy++) {
  84     for (ix = 0; ix < nx; ix++)
  85       pRow[ix] = vReal[ix][iy];
  86     Fourier::shuffleNaturalToFourierOrder (pRow, nx);
  87     for (ix = 0; ix < nx; ix++)
  88       vReal[ix][iy] = pRow[ix];
  89     if (im.isComplex()) {
  90       for (ix = 0; ix < nx; ix++)
  91         pRow[ix] = vImag[ix][iy];
  92       Fourier::shuffleNaturalToFourierOrder (pRow, nx);
  93       for (ix = 0; ix < nx; ix++)
  94         vImag[ix][iy] = pRow[ix];
  95     }
  96   }
  97   delete [] pRow;
  98 }
  99
 100 #ifdef HAVE_FFTW
 101 void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
 102 {
 103   fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
 104   int i;
 105
 106   if (isOdd(n)) { // Odd
 107     int iHalfN = (n - 1) / 2;
 108
 109     pTemp[0][0] = pVector[iHalfN][0];
 110     pTemp[0][1] = pVector[iHalfN][1];
 111     for (i = 0; i < iHalfN; i++) {
 112       pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
 113       pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
 114     }
 115     for (i = 0; i < iHalfN; i++) {
 116       pTemp[i + iHalfN + 1][0] = pVector[i][0];
 117       pTemp[i + iHalfN + 1][1] = pVector[i][1];
 118     }
 119   } else {     // Even
 120     int iHalfN = n / 2;
 121     pTemp[0][0] = pVector[iHalfN][0];
 122     pTemp[0][1] = pVector[iHalfN][1];
 123     for (i = 0; i < iHalfN - 1; i++) {
 124       pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
 125       pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
 126     }
 127     for (i = 0; i < iHalfN; i++) {
 128       pTemp[i + iHalfN][0] = pVector[i][0];
 129       pTemp[i + iHalfN][1] = pVector[i][1];
 130     }
 131   }
 132
 133   for (i = 0; i < n; i++) {
 134     pVector[i][0] = pTemp[i][0];
 135     pVector[i][1] = pTemp[i][1];
 136   }
 137   fftw_free(pTemp);
 138 }
 139
 140 void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
 141 {
 142   fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
 143   int i;
 144   if (isOdd(n)) { // Odd
 145     int iHalfN = (n - 1) / 2;
 146
 147     pTemp[iHalfN][0] = pVector[0][0];
 148     pTemp[iHalfN][1] = pVector[0][1];
 149     for (i = 0; i < iHalfN; i++) {
 150       pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
 151       pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
 152     }
 153     for (i = 0; i < iHalfN; i++) {
 154       pTemp[i][0] = pVector[i + iHalfN + 1][0];
 155       pTemp[i][1] = pVector[i + iHalfN + 1][1];
 156     }
 157   } else {     // Even
 158     int iHalfN = n / 2;
 159     pTemp[iHalfN][0] = pVector[0][0];
 160     pTemp[iHalfN][1] = pVector[0][1];
 161     for (i = 0; i < iHalfN; i++) {
 162       pTemp[i][0] = pVector[i + iHalfN][0];
 163       pTemp[i][1] = pVector[i + iHalfN][1];
 164     }
 165     for (i = 0; i < iHalfN - 1; i++) {
 166       pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
 167       pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
 168     }
 169   }
 170
 171   for (i = 0; i < n; i++) {
 172     pVector[i][0] = pTemp[i][0];
 173     pVector[i][1] = pTemp[i][1];
 174   }
 175
 176   fftw_free(pTemp);
 177 }
 178 #endif
 179