X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffourier.cpp;h=d85bfbe5b97d8c332a582ed94eea94469aaca171;hp=06145a3c628495a3736ea7c3bda724beedff2020;hb=f13a8c004b8f182b42d9e4df2bcd7c7f030bf1ad;hpb=9f29c8b32c972db1178d6f8551d5cd57ceb67083 diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp index 06145a3..d85bfbe 100644 --- a/libctsim/fourier.cpp +++ b/libctsim/fourier.cpp @@ -7,9 +7,7 @@ ** Date Started: Dec 2000 ** ** 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 $ +** Copyright (c) 1983-2009 Kevin Rosenberg ** ** 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 +97,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