Merge pull request #160 from wjp/FBP_filters

Custom filter support for CPU FBP

2 files changed, 28 insertions, 423 deletions

diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index 48fb7dc..a5b8a7a 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -35,27 +35,18 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/fdk.h"
 
 #include "astra/Logging.h"
+#include "astra/Filters.h"
 
 #include <cuda.h>
 
 namespace astraCUDA {
 
 
-
-static int calcNextPowerOfTwo(int n)
-{
-	int x = 1;
-	while (x < n)
-		x *= 2;
-
-	return x;
-}
-
 // static
 int FBP::calcFourierFilterSize(int _iDetectorCount)
 {
-	int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount);
-	int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+	int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * _iDetectorCount);
+	int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 	// CHECKME: Matlab makes this at least 64. Do we also need to?
 	return iFreqBinCount;
@@ -88,7 +79,7 @@ bool FBP::init()
 	return true;
 }
 
-bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */)
+bool FBP::setFilter(const astra::SFilterConfig &_cfg)
 {
 	if (D_filter)
 	{
@@ -96,19 +87,19 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 		D_filter = 0;
 	}
 
-	if (_eFilter == astra::FILTER_NONE)
+	if (_cfg.m_eType == astra::FILTER_NONE)
 		return true; // leave D_filter set to 0
 
 
-	int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
-	int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+	int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
+	int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 	cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount];
 	memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount);
 
 	allocateComplexOnDevice(dims.iProjAngles, iFreqBinCount, (cufftComplex**)&D_filter);
 
-	switch(_eFilter)
+	switch(_cfg.m_eType)
 	{
 		case astra::FILTER_NONE:
 			// handled above
@@ -130,7 +121,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 		case astra::FILTER_FLATTOP:
 		case astra::FILTER_PARZEN:
 		{
-			genFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter);
+			genCuFFTFilter(_cfg, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount);
 			uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter);
 
 			break;
@@ -138,11 +129,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 		case astra::FILTER_PROJECTION:
 		{
 			// make sure the offered filter has the correct size
-			assert(_iFilterWidth == iFreqBinCount);
+			assert(_cfg.m_iCustomFilterWidth == iFreqBinCount);
 
 			for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++)
 			{
-				float fValue = _pfHostFilter[iFreqBinIndex];
+				float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex];
 
 				for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
 				{
@@ -156,13 +147,13 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 		case astra::FILTER_SINOGRAM:
 		{
 			// make sure the offered filter has the correct size
-			assert(_iFilterWidth == iFreqBinCount);
+			assert(_cfg.m_iCustomFilterWidth == iFreqBinCount);
 
 			for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++)
 			{
 				for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
 				{
-					float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth];
+					float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth];
 
 					pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue;
 					pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f;
@@ -178,16 +169,16 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 			float * pfHostRealFilter = new float[iRealFilterElementCount];
 			memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount);
 
-			int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount);
-			int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2;
+			int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount);
+			int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
 			int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
 
-			int iFilterShiftSize = _iFilterWidth / 2;
+			int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2;
 
 			for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++)
 			{
 				int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount;
-				float fValue = _pfHostFilter[iDetectorIndex];
+				float fValue = _cfg.m_pfCustomFilter[iDetectorIndex];
 
 				for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
 				{
@@ -213,11 +204,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 			float* pfHostRealFilter = new float[iRealFilterElementCount];
 			memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount);
 
-			int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount);
-			int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2;
+			int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount);
+			int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
 			int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
 
-			int iFilterShiftSize = _iFilterWidth / 2;
+			int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2;
 
 			for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++)
 			{
@@ -225,7 +216,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
 
 				for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
 				{
-					float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth];
+					float fValue = _cfg.m_pfCustomFilter[iDetectorIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth];
 					pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue;
 				}
 			}
@@ -310,8 +301,8 @@ bool FBP::iterate(unsigned int iterations)
 
 	if (D_filter) {
 
-		int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
-		int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount);
+		int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
+		int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 		cufftComplex * pDevComplexSinogram = NULL;
 
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index bd8cab5..2e94b79 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -275,17 +275,6 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex,
 	return true;
 }
 
-
-// Because the input is real, the Fourier transform is symmetric.
-// CUFFT only outputs the first half (ignoring the redundant second half),
-// and expects the same as input for the IFFT.
-int calcFFTFourierSize(int _iFFTRealSize)
-{
-	int iFFTFourierSize = _iFFTRealSize / 2 + 1;
-
-	return iFFTFourierSize;
-}
-
 void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter,
                    int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount)
 {
@@ -300,387 +289,13 @@ void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter,
 	}
 }
 
-void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
+void genCuFFTFilter(const SFilterConfig &_cfg, int _iProjectionCount,
                cufftComplex * _pFilter, int _iFFTRealDetectorCount,
-               int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */)
+               int _iFFTFourierDetectorCount)
 {
-	float * pfFilt = new float[_iFFTFourierDetectorCount];
-	float * pfW = new float[_iFFTFourierDetectorCount];
-
-	// We cache one Fourier transform for repeated FBP's of the same size
-	static float *pfData = 0;
-	static int iFilterCacheSize = 0;
-
-	if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) {
-		// Compute filter in spatial domain
-
-		delete[] pfData;
-		pfData = new float[2*_iFFTRealDetectorCount];
-		int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)];
-		ip[0] = 0;
-		float32 *w = new float32[_iFFTRealDetectorCount/2];
-
-		for (int i = 0; i < _iFFTRealDetectorCount; ++i) {
-			pfData[2*i+1] = 0.0f;
-
-			if (i & 1) {
-				int j = i;
-				if (2*j > _iFFTRealDetectorCount)
-					j = _iFFTRealDetectorCount - j;
-				float f = M_PI * j;
-				pfData[2*i] = -1 / (f*f);
-			} else {
-				pfData[2*i] = 0.0f;
-			}
-		}
-
-		pfData[0] = 0.25f;
-
-		cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
-		delete[] ip;
-		delete[] w;
-
-		iFilterCacheSize = _iFFTRealDetectorCount;
-	}
-
-	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-	{
-		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
-
-		pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex];
-		pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
-	}
-
-	switch(_eFilter)
-	{
-		case FILTER_RAMLAK:
-		{
-			// do nothing
-			break;
-		}
-		case FILTER_SHEPPLOGAN:
-		{
-			// filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)))
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD));
-			}
-			break;
-		}
-		case FILTER_COSINE:
-		{
-			// filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d))
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _fD);
-			}
-			break;
-		}
-		case FILTER_HAMMING:
-		{
-			// filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d))
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _fD));
-			}
-			break;
-		}
-		case FILTER_HANN:
-		{
-			// filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _fD)) / 2.0f;
-			}
-			break;
-		}
-		case FILTER_TUKEY:
-		{
-			float fAlpha = _fParameter;
-			if(_fParameter < 0.0f) fAlpha = 0.5f;
-			float fN = (float)_iFFTFourierDetectorCount;
-			float fHalfN = fN / 2.0f;
-			float fEnumTerm = fAlpha * fHalfN;
-			float fDenom = (1.0f - fAlpha) * fHalfN;
-			float fBlockStart = fHalfN - fEnumTerm;
-			float fBlockEnd = fHalfN + fEnumTerm;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fAbsSmallN = fabs((float)iDetectorIndex);
-				float fStoredValue = 0.0f;
-
-				if((fBlockStart <= fAbsSmallN) && (fAbsSmallN <= fBlockEnd))
-				{
-					fStoredValue = 1.0f;
-				}
-				else
-				{
-					float fEnum = fAbsSmallN - fEnumTerm;
-					float fCosInput = M_PI * fEnum / fDenom;
-					fStoredValue = 0.5f * (1.0f + cosf(fCosInput));
-				}
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_LANCZOS:
-		{
-			float fDenum = (float)(_iFFTFourierDetectorCount - 1);
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fX = 2.0f * fSmallN / fDenum - 1.0f;
-				float fSinInput = M_PI * fX;
-				float fStoredValue = 0.0f;
-
-				if(fabsf(fSinInput) > 0.001f)
-				{
-					fStoredValue = sin(fSinInput)/fSinInput;
-				}
-				else
-				{
-					fStoredValue = 1.0f;
-				}
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_TRIANGULAR:
-		{
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1);
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fAbsInput = fSmallN - fNMinusOne / 2.0f;
-				float fParenInput = fNMinusOne / 2.0f - fabsf(fAbsInput);
-				float fStoredValue = 2.0f / fNMinusOne * fParenInput;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_GAUSSIAN:
-		{
-			float fSigma = _fParameter;
-			if(_fParameter < 0.0f) fSigma = 0.4f;
-			float fN = (float)_iFFTFourierDetectorCount;
-			float fQuotient = (fN - 1.0f) / 2.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fEnum = fSmallN - fQuotient;
-				float fDenom = fSigma * fQuotient;
-				float fPower = -0.5f * (fEnum / fDenom) * (fEnum / fDenom);
-				float fStoredValue = expf(fPower);
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_BARTLETTHANN:
-		{
-			const float fA0 = 0.62f;
-			const float fA1 = 0.48f;
-			const float fA2 = 0.38f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fAbsInput = fSmallN / fNMinusOne - 0.5f;
-				float fFirstTerm = fA1 * fabsf(fAbsInput);
-				float fCosInput = 2.0f * M_PI * fSmallN / fNMinusOne;
-				float fSecondTerm = fA2 * cosf(fCosInput);
-				float fStoredValue = fA0 - fFirstTerm - fSecondTerm;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_BLACKMAN:
-		{
-			float fAlpha = _fParameter;
-			if(_fParameter < 0.0f) fAlpha = 0.16f;
-			float fA0 = (1.0f - fAlpha) / 2.0f;
-			float fA1 = 0.5f;
-			float fA2 = fAlpha / 2.0f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1);
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fCosInput1 = 2.0f * M_PI * 0.5f * fSmallN / fNMinusOne;
-				float fCosInput2 = 4.0f * M_PI * 0.5f * fSmallN / fNMinusOne;
-				float fStoredValue = fA0 - fA1 * cosf(fCosInput1) + fA2 * cosf(fCosInput2);
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_NUTTALL:
-		{
-			const float fA0 = 0.355768f;
-			const float fA1 = 0.487396f;
-			const float fA2 = 0.144232f;
-			const float fA3 = 0.012604f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
-				float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
-				float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
-				float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
-				float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_BLACKMANHARRIS:
-		{
-			const float fA0 = 0.35875f;
-			const float fA1 = 0.48829f;
-			const float fA2 = 0.14128f;
-			const float fA3 = 0.01168f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
-				float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
-				float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
-				float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
-				float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_BLACKMANNUTTALL:
-		{
-			const float fA0 = 0.3635819f;
-			const float fA1 = 0.4891775f;
-			const float fA2 = 0.1365995f;
-			const float fA3 = 0.0106411f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
-				float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
-				float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
-				float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
-				float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_FLATTOP:
-		{
-			const float fA0 = 1.0f;
-			const float fA1 = 1.93f;
-			const float fA2 = 1.29f;
-			const float fA3 = 0.388f;
-			const float fA4 = 0.032f;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
-				float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
-				float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
-				float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
-				float fFourthTerm = fA4 * cosf(8.0f * fBaseCosInput);
-				float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm + fFourthTerm;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_KAISER:
-		{
-			float fAlpha = _fParameter;
-			if(_fParameter < 0.0f) fAlpha = 3.0f;
-			float fPiTimesAlpha = M_PI * fAlpha;
-			float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1);
-			float fDenom = (float)j0((double)fPiTimesAlpha);
-
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fSquareInput = 2.0f * fSmallN / fNMinusOne - 1;
-				float fSqrtInput = 1.0f - fSquareInput * fSquareInput;
-				float fBesselInput = fPiTimesAlpha * sqrt(fSqrtInput);
-				float fEnum = (float)j0((double)fBesselInput);
-				float fStoredValue = fEnum / fDenom;
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		case FILTER_PARZEN:
-		{
-			for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-			{
-				float fSmallN = (float)iDetectorIndex;
-				float fQ = fSmallN / (float)(_iFFTFourierDetectorCount - 1);
-				float fStoredValue = 0.0f;
-
-				if(fQ <= 0.5f)
-				{
-					fStoredValue = 1.0f - 6.0f * fQ * fQ * (1.0f - fQ);
-				}
-				else
-				{
-					float fCubedValue = 1.0f - fQ;
-					fStoredValue = 2.0f * fCubedValue * fCubedValue * fCubedValue;
-				}
-
-				pfFilt[iDetectorIndex] *= fStoredValue;
-			}
-
-			break;
-		}
-		default:
-		{
-			ASTRA_ERROR("Cannot serve requested filter");
-		}
-	}
-
-	// filt(w>pi*d) = 0;
-	float fPiTimesD = M_PI * _fD;
-	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-	{
-		float fWValue = pfW[iDetectorIndex];
-
-		if(fWValue > fPiTimesD)
-		{
-			pfFilt[iDetectorIndex] = 0.0f;
-		}
-	}
+	float * pfFilt = astra::genFilter(_cfg,
+	                                  _iFFTRealDetectorCount,
+	                                  _iFFTFourierDetectorCount);
 
 	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
 	{
@@ -695,7 +310,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 	}
 
 	delete[] pfFilt;
-	delete[] pfW;
 }