X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffourier.cpp;h=08940d43636881e2db2d632b09ccfba95fae1ac3;hp=09033acd3cecb1ad1478738f244256f0194b2353;hb=c00c639073653fac7463a88f2b000f263236550d;hpb=23b7ef994fc5d95fcca6d4ae69abbd5971101262 diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp index 09033ac..08940d4 100644 --- a/libctsim/fourier.cpp +++ b/libctsim/fourier.cpp @@ -1,279 +1,279 @@ -/***************************************************************************** -** FILE IDENTIFICATION -** -** Name: fourier.cpp -** Purpose: Fourier transform functions -** Programmer: Kevin Rosenberg -** Date Started: Dec 2000 -** -** This is part of the CTSim program -** Copyright (C) 1983-2001 Kevin Rosenberg -** -** $Id: fourier.cpp,v 1.1 2001/01/02 06:33:04 kevin Exp $ -** -** This program is free software; you can redistribute it and/or modify -** it under the terms of the GNU General Public License (version 2) as -** published by the Free Software Foundation. -** -** This program is distributed in the hope that it will be useful, -** but WITHOUT ANY WARRANTY; without even the implied warranty of -** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -** GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License -** along with this program; if not, write to the Free Software -** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -******************************************************************************/ - -#include "ct.h" - - - -void -Fourier::shuffleFourierToNaturalOrder (ImageFile& im) -{ - ImageFileArray vReal = im.getArray(); - ImageFileArray vImag = im.getImaginaryArray(); - unsigned int ix, iy; - unsigned int nx = im.nx(); - unsigned int ny = im.ny(); - - // shuffle each column - for (ix = 0; ix < nx; ix++) { - Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny); - if (im.isComplex()) - Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny); - } - - // shuffle each row - float* pRow = new float [nx]; - for (iy = 0; iy < ny; iy++) { - for (ix = 0; ix < nx; ix++) - pRow[ix] = vReal[ix][iy]; - Fourier::shuffleFourierToNaturalOrder (pRow, nx); - for (ix = 0; ix < nx; ix++) - vReal[ix][iy] = pRow[ix]; - if (im.isComplex()) { - for (ix = 0; ix < nx; ix++) - pRow[ix] = vImag[ix][iy]; - Fourier::shuffleFourierToNaturalOrder (pRow, nx);; - for (ix = 0; ix < nx; ix++) - vImag[ix][iy] = pRow[ix]; - } - } - delete pRow; -} - -void -Fourier::shuffleNaturalToFourierOrder (ImageFile& im) -{ - ImageFileArray vReal = im.getArray(); - ImageFileArray vImag = im.getImaginaryArray(); - unsigned int ix, iy; - unsigned int nx = im.nx(); - unsigned int ny = im.ny(); - - // shuffle each x column - for (ix = 0; ix < nx; ix++) { - Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny); - if (im.isComplex()) - Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny); - } - - // shuffle each y row - float* pRow = new float [nx]; - for (iy = 0; iy < ny; iy++) { - for (ix = 0; ix < nx; ix++) - pRow[ix] = vReal[ix][iy]; - Fourier::shuffleNaturalToFourierOrder (pRow, nx); - for (ix = 0; ix < nx; ix++) - vReal[ix][iy] = pRow[ix]; - if (im.isComplex()) { - for (ix = 0; ix < nx; ix++) - pRow[ix] = vImag[ix][iy]; - Fourier::shuffleNaturalToFourierOrder (pRow, nx); - for (ix = 0; ix < nx; ix++) - vImag[ix][iy] = pRow[ix]; - } - } - delete [] pRow; -} - - -// Odd Number of Points -// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2 -// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1 -// Even Number of Points -// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) -// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 - -void -Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n) -{ - double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + 1] = pdVector[i + iHalfN + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN] = pdVector[i]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - -void -Fourier::shuffleNaturalToFourierOrder (std::complex* pdVector, const int n) -{ - std::complex* pdTemp = new std::complex [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + 1] = pdVector[i + iHalfN + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN] = pdVector[i]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete [] pdTemp; -} - - -void -Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n) -{ - float* pdTemp = new float [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + 1] = pdVector[i + iHalfN + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN] = pdVector[i]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - - - -void -Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n) -{ - double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i+1]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - - -void -Fourier::shuffleFourierToNaturalOrder (std::complex* pdVector, const int n) -{ - std::complex* pdTemp = new std::complex [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i+1]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete [] pdTemp; -} - - - - -void -Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n) -{ - float* pTemp = new float [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pTemp[iHalfN] = pVector[0]; - for (i = 0; i < iHalfN; i++) - pTemp[i + 1 + iHalfN] = pVector[i + 1]; - for (i = 0; i < iHalfN; i++) - pTemp[i] = pVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pTemp[iHalfN] = pVector[0]; - for (i = 0; i < iHalfN; i++) - pTemp[i] = pVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pTemp[i + iHalfN + 1] = pVector[i+1]; - } - - for (i = 0; i < n; i++) - pVector[i] = pTemp[i]; - delete [] pTemp; -} - - +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: fourier.cpp +** Purpose: Fourier transform functions +** Programmer: Kevin Rosenberg +** Date Started: Dec 2000 +** +** This is part of the CTSim program +** Copyright (C) 1983-2001 Kevin Rosenberg +** +** $Id: fourier.cpp,v 1.2 2001/01/02 16:02:13 kevin Exp $ +** +** This program is free software; you can redistribute it and/or modify +** it under the terms of the GNU General Public License (version 2) as +** published by the Free Software Foundation. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program; if not, write to the Free Software +** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +******************************************************************************/ + +#include "ct.h" + + + +void +Fourier::shuffleFourierToNaturalOrder (ImageFile& im) +{ + ImageFileArray vReal = im.getArray(); + ImageFileArray vImag = im.getImaginaryArray(); + unsigned int ix, iy; + unsigned int nx = im.nx(); + unsigned int ny = im.ny(); + + // shuffle each column + for (ix = 0; ix < nx; ix++) { + Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny); + if (im.isComplex()) + Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny); + } + + // shuffle each row + float* pRow = new float [nx]; + for (iy = 0; iy < ny; iy++) { + for (ix = 0; ix < nx; ix++) + pRow[ix] = vReal[ix][iy]; + Fourier::shuffleFourierToNaturalOrder (pRow, nx); + for (ix = 0; ix < nx; ix++) + vReal[ix][iy] = pRow[ix]; + if (im.isComplex()) { + for (ix = 0; ix < nx; ix++) + pRow[ix] = vImag[ix][iy]; + Fourier::shuffleFourierToNaturalOrder (pRow, nx);; + for (ix = 0; ix < nx; ix++) + vImag[ix][iy] = pRow[ix]; + } + } + delete pRow; +} + +void +Fourier::shuffleNaturalToFourierOrder (ImageFile& im) +{ + ImageFileArray vReal = im.getArray(); + ImageFileArray vImag = im.getImaginaryArray(); + unsigned int ix, iy; + unsigned int nx = im.nx(); + unsigned int ny = im.ny(); + + // shuffle each x column + for (ix = 0; ix < nx; ix++) { + Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny); + if (im.isComplex()) + Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny); + } + + // shuffle each y row + float* pRow = new float [nx]; + for (iy = 0; iy < ny; iy++) { + for (ix = 0; ix < nx; ix++) + pRow[ix] = vReal[ix][iy]; + Fourier::shuffleNaturalToFourierOrder (pRow, nx); + for (ix = 0; ix < nx; ix++) + vReal[ix][iy] = pRow[ix]; + if (im.isComplex()) { + for (ix = 0; ix < nx; ix++) + pRow[ix] = vImag[ix][iy]; + Fourier::shuffleNaturalToFourierOrder (pRow, nx); + for (ix = 0; ix < nx; ix++) + vImag[ix][iy] = pRow[ix]; + } + } + delete [] pRow; +} + + +// Odd Number of Points +// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2 +// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1 +// Even Number of Points +// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) +// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 + +void +Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n) +{ + double* pdTemp = new double [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + 1] = pdVector[i + iHalfN + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + +void +Fourier::shuffleNaturalToFourierOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + 1] = pdVector[i + iHalfN + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} + + +void +Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n) +{ + float* pdTemp = new float [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + 1] = pdVector[i + iHalfN + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + + + +void +Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n) +{ + double* pdTemp = new double [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + + +void +Fourier::shuffleFourierToNaturalOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} + + + + +void +Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n) +{ + float* pTemp = new float [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pTemp[iHalfN] = pVector[0]; + for (i = 0; i < iHalfN; i++) + pTemp[i + 1 + iHalfN] = pVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pTemp[i] = pVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pTemp[iHalfN] = pVector[0]; + for (i = 0; i < iHalfN; i++) + pTemp[i] = pVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pTemp[i + iHalfN + 1] = pVector[i+1]; + } + + for (i = 0; i < n; i++) + pVector[i] = pTemp[i]; + delete [] pTemp; +} + +