From b2fc6c70434674d74551c3a6c01ffb3233499312 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 1 Jul 2013 22:34:11 +0000 Subject: Update version to 1.3 --- src/Fourier.cpp | 233 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 233 insertions(+) create mode 100644 src/Fourier.cpp (limited to 'src/Fourier.cpp') diff --git a/src/Fourier.cpp b/src/Fourier.cpp new file mode 100644 index 0000000..0f7da28 --- /dev/null +++ b/src/Fourier.cpp @@ -0,0 +1,233 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +#include "astra/Fourier.h" + +namespace astra { + + +void discreteFourierTransform1D(unsigned int iLength, + const float32* pfRealIn, + const float32* pfImaginaryIn, + float32* pfRealOut, + float32* pfImaginaryOut, + unsigned int iStrideIn, + unsigned int iStrideOut, + bool inverse) +{ + for (unsigned int w = 0; w < iLength; w++) + { + pfRealOut[iStrideOut*w] = pfImaginaryOut[iStrideOut*w] = 0; + for (unsigned int y = 0; y < iLength; y++) + { + float32 a = 2 * PI * w * y / float32(iLength); + if (!inverse) + a = -a; + float32 ca = cos(a); + float32 sa = sin(a); + pfRealOut[iStrideOut*w] += pfRealIn[iStrideIn*y] * ca - pfImaginaryIn[iStrideIn*y] * sa; + pfImaginaryOut[iStrideOut*w] += pfRealIn[iStrideIn*y] * sa + pfImaginaryIn[iStrideIn*y] * ca; + } + } + + if (inverse) { + for (unsigned int x = 0; x < iLength; ++x) { + pfRealOut[iStrideOut*x] /= iLength; + pfImaginaryOut[iStrideOut*x] /= iLength; + } + } +} + +void discreteFourierTransform2D(unsigned int iHeight, unsigned int iWidth, + const float32* pfRealIn, + const float32* pfImaginaryIn, + float32* pfRealOut, + float32* pfImaginaryOut, + bool inverse) +{ + float32* reTemp = new float32[iWidth * iHeight]; + float32* imTemp = new float32[iWidth * iHeight]; + + //calculate the fourier transform of the columns + for (unsigned int x = 0; x < iWidth; x++) + { + discreteFourierTransform1D(iHeight, pfRealIn+x, pfImaginaryIn+x, + reTemp+x, imTemp+x, + iWidth, iWidth, inverse); + } + + //calculate the fourier transform of the rows + for(unsigned int y = 0; y < iHeight; y++) + { + discreteFourierTransform1D(iWidth, + reTemp+y*iWidth, + imTemp+y*iWidth, + pfRealOut+y*iWidth, + pfImaginaryOut+y*iWidth, + 1, 1, inverse); + } + + delete[] reTemp; + delete[] imTemp; +} + +/** permute the entries from pfDataIn into pfDataOut to prepare for an + * in-place FFT. pfDataIn may be equal to pfDataOut. + */ +static void bitReverse(unsigned int iLength, + const float32* pfDataIn, float32* pfDataOut, + unsigned int iStrideShiftIn, + unsigned int iStrideShiftOut) +{ + if (pfDataIn == pfDataOut) { + assert(iStrideShiftIn == iStrideShiftOut); + float32 t; + unsigned int j = 0; + for(unsigned int i = 0; i < iLength - 1; i++) { + if (i < j) { + t = pfDataOut[i< 1) { + n /= 2; + ++l; + } + return l; +} + +/** perform 1D FFT. iLength, iStrideIn, iStrideOut must be powers of two. */ +void fastTwoPowerFourierTransform1D(unsigned int iLength, + const float32* pfRealIn, + const float32* pfImaginaryIn, + float32* pfRealOut, + float32* pfImaginaryOut, + unsigned int iStrideIn, + unsigned int iStrideOut, + bool inverse) +{ + unsigned int iStrideShiftIn = log2(iStrideIn); + unsigned int iStrideShiftOut = log2(iStrideOut); + unsigned int iLogLength = log2(iLength); + + bitReverse(iLength, pfRealIn, pfRealOut, iStrideShiftIn, iStrideShiftOut); + bitReverse(iLength, pfImaginaryIn, pfImaginaryOut, iStrideShiftIn, iStrideShiftOut); + + float32 ca = -1.0; + float32 sa = 0.0; + unsigned int l1 = 1, l2 = 1; + for(unsigned int l=0; l < iLogLength; ++l) + { + l1 = l2; + l2 *= 2; + float32 u1 = 1.0; + float32 u2 = 0.0; + for(unsigned int j = 0; j < l1; j++) + { + for(unsigned int i = j; i < iLength; i += l2) + { + unsigned int i1 = i + l1; + float32 t1 = u1 * pfRealOut[i1<