X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffourier.cpp;h=149b4de234b5fc2c4738295ce04df5c51e792dd6;hp=06145a3c628495a3736ea7c3bda724beedff2020;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=9f29c8b32c972db1178d6f8551d5cd57ceb67083 diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp index 06145a3..149b4de 100644 --- a/libctsim/fourier.cpp +++ b/libctsim/fourier.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: fourier.cpp,v 1.4 2001/01/28 19:10:18 kevin Exp $ +** $Id$ ** ** 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 @@ -99,181 +99,83 @@ Fourier::shuffleNaturalToFourierOrder (ImageFile& im) 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) +#ifdef HAVE_FFTW +void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n) { - double* pdTemp = new double [n]; + fftw_complex* pTemp = static_cast(fftw_malloc(sizeof(fftw_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 (std::complex* pdVector, const int n) -{ - std::complex* pdTemp = new std::complex [n]; - int i; - if (n % 2) { // Odd + if (isOdd(n)) { // 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]; + pTemp[0][0] = pVector[iHalfN][0]; + pTemp[0][1] = pVector[iHalfN][1]; + for (i = 0; i < iHalfN; i++) { + pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0]; + pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1]; + } + for (i = 0; i < iHalfN; i++) { + pTemp[i + iHalfN + 1][0] = pVector[i][0]; + pTemp[i + iHalfN + 1][1] = pVector[i][1]; + } } 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]; + pTemp[0][0] = pVector[iHalfN][0]; + pTemp[0][1] = pVector[iHalfN][1]; + for (i = 0; i < iHalfN - 1; i++) { + pTemp[i + 1][0] = pVector[i + iHalfN + 1][0]; + pTemp[i + 1][1] = pVector[i + iHalfN + 1][1]; + } + for (i = 0; i < iHalfN; i++) { + pTemp[i + iHalfN][0] = pVector[i][0]; + pTemp[i + iHalfN][1] = pVector[i][1]; + } } - 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++) { + pVector[i][0] = pTemp[i][0]; + pVector[i][1] = pTemp[i][1]; } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; + fftw_free(pTemp); } - - -void -Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n) +void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n) { - double* pdTemp = new double [n]; + fftw_complex* pTemp = static_cast(fftw_malloc(sizeof(fftw_complex) * n)); int i; - if (n % 2) { // Odd + if (isOdd(n)) { // 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]; + pTemp[iHalfN][0] = pVector[0][0]; + pTemp[iHalfN][1] = pVector[0][1]; + for (i = 0; i < iHalfN; i++) { + pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0]; + pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1]; + } + for (i = 0; i < iHalfN; i++) { + pTemp[i][0] = pVector[i + iHalfN + 1][0]; + pTemp[i][1] = pVector[i + iHalfN + 1][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]; + pTemp[iHalfN][0] = pVector[0][0]; + pTemp[iHalfN][1] = pVector[0][1]; + for (i = 0; i < iHalfN; i++) { + pTemp[i][0] = pVector[i + iHalfN][0]; + pTemp[i][1] = pVector[i + iHalfN][1]; + } + for (i = 0; i < iHalfN - 1; i++) { + pTemp[i + iHalfN + 1][0] = pVector[i+1][0]; + pTemp[i + iHalfN + 1][1] = pVector[i+1][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][0] = pTemp[i][0]; + pVector[i][1] = pTemp[i][1]; } - - for (i = 0; i < n; i++) - pVector[i] = pTemp[i]; - delete [] pTemp; -} + fftw_free(pTemp); +} +#endif